Quantifying Fluctuation/Correlation Effects on the Order-Disorder Transition of 

Symmetric Diblock Copolymers 
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Using fast off-lattice Monte Carlo simulations with experimentally accessible fluctuations, we re- 
port the first systematic study unambiguously quantifying the shift of the order-disorder transition 
(ODT) x* °f symmetric diblock copolymers from the mean-field prediction Xmf- 0m simulations 
are performed in a canonical ensemble with variable box lengths to eliminate the restriction of pe- 
riodic boundary conditions on the lamellar period. Exactly the same model system (Hamiltonian) 
is used in both our simulations and mean-field theory; the ODT shift is therefore due to the fluc- 
tuations/correlations neglected by the latter. While X*/Xmf — 1 oc N~ k is found with tf denoting 
the invariant degree of polymerization, k decreases around the A/"-value corresponding to the close 
packing of polymer segments as hard spheres, indicating the short-range correlation effects. 

PACS numbers: 64.60.Cn, 64.70.km 



Owing to the well developed polymer self-consistent 
field (SCF) calculations 0, good or even quantitative un- 
derstanding has been achieved for the self-assembly of 
flexible linear diblock copolymer (DBC) melts in bulkQ. 
Due to its mean-field approximation, however, SCF the- 
ory gives qualitatively incorrect predictions in the region 
near the order-disorder transition (ODT) where the sys- 
tem fluctuations it neglects become important. In par- 
ticular, it fails to capture the fluctuation-induced first- 
order phase transition for ODT of symmetric DBC0 and 
the direct transition between the gyroid and disordered 
phases Q. Although the former is a classic problem in 
polymer science and has been extensively studied by ex- 
periments, theories and simulations, no quantitative un- 
derstanding of the fluctuation effects on ODT is achieved. 
Fig. [T] summarizes recent simulation results on the ODT 
shift of symmetric DBC from the mean-field prediction 
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FIG. 1: (Color online) Some literature results on the ODT 
shift. See text for details. 
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vs. the invariant degree of polymerization M = nB? e q/V, 
where n denotes the number of copolymer chains, i? e o 
the end-to-end distance of an ideal chain, and V the sys- 
tem volume; these are the most accurate data obtained 
using each method as explained below. 

Beardsley and Matsen performed conventional Monte 
Carlo (MC) simulations on a face-centered cubic (FCC) 
lattice in a canonical ensemble with replica-exchange to 
study symmetric DBC of N — 30 segments modeled 
by the self- and mutual-avoiding walks with nearest- 
neighbor repulsion between A and B segments cab (in 
units of ksT, where ks is the Boltzmann constant and 
T the thermodynamic temperature); about 20% lattice 
sites were unoccupied, treated as an athermal solvent. @ 
From the peak of constant-volume heat capacity, they 
determined ODT to be \* N = ze* AB N = 40.5, where 
z = 12 is the lattice coordination number. [j| The mean- 
field ODT was determined to be Xmf-^ = 14-654 using 
lattice SCF calculations based on the same Hamiltonian 
as used in their simulations. Q 

Using fast off-lattice Monte Carlo (FOMC) simula- 
tions in an isothermal-isobaric ensemble, de Pablo and 
co-workers studied compressible symmetric DBC melts 
modeled by discrete Gaussian chains (DGC) of N = 64 
with a position-independent but anisotropic pair poten- 
tial of cubic symmetry. Q They determined x*N by cal- 
culating the chemical potentials of the disordered and 
lamellar phases. Q Note that their x*A was renormalized 
to take into account the short-range correlations (due to 
the finite interaction range used in their simulations) ap- 
proximately; their Xmf-^ was therefore determined using 
the random-phase approximation for DGC with Dirac 5- 
function interactions. 0, H[ 

Muller and Daoulas performed single-chain-in-mean- 
field (SCMF) simulations to study compressible symmet- 
ric DBC melts modeled by DGC of N = 32, and deter- 
mined x* A = 13.65 ±0.1 by comparing the free energies 
of the disordered and lamellar phases. @ Their SCMF 
simulation is similar to canonical-ensemble FOMC sim- 
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ulation with a spatial discretization scheme (i.e., with 
a position-dependent and anisotropic pair potential) [To| , 
except that a second-order term in the energy differ- 
ence due to MC trial moves is neglected by the quasi- 
instantaneous field approximation[llj. They also esti- 
mated Xmf^ — 11.7 ± 0.5 by fitting the lamellar com- 
position profiles (averaged over directions parallel to the 
lamellar interfaces) of a system at Af w 2.5 x 10 9 at var- 
ious xiV to a sinusoidal function and then extrapolating 
the square of the so-obtained amplitude of the composi- 
tion profiles to 0. 9] 

Using field-theoretic simulations in a canonical en- 
semble, Fredrickson and co-workers studied incompress- 
ible DBC melts modeled by continuous Gaussian chains 
with Dirac 6- function interactions. [l2| They determined 
X*N = 12.2 (for slightly asymmetric DBC with a com- 
position / = 0.49) by comparing the free energies of the 
disordered and lamellar phases. [12j For this "standard" 
model (which cannot be directly used/tested in molec- 
ular simulations), Xmf^ = 10.495 for / = 0.5 is the 
well-known Leibler's result [l3[ and is used in Fig. [TJ 

Finally, the prediction of the fluctuation theory[14| by 
Fredrickson and Helfand (FH), x*/xm F -1 = 3.91A/"- 1 / 3 , 
is also shown as a straight line in Fig. [1] This theory is 
based on the Ohta-Kawasaki effective Hamiltonian [15| for 
the above "standard" model and the Hartree analysis by 
Brazovskii|l6|. Although the latter is rigorously accurate 



only for Af > 10 10 [l4|, the FH prediction has often been 



compared with simulation results at much smaller Af in 
the literature. With the different models and methods 
used and the scarce simulation data in Fig. [TJ it is clear 
that the ODT shift of symmetric DBC due to fluctuations 
is a problem far from being well understood. 

The system fluctuations are controlled by AT, as sug- 
gested by the FH theory. For monodisperse DBC melts 
where each copolymer chain consists of N A monomers 
of type A followed by Nb monomers of type B, with 



N A a 2 A + N B a B and nN — p V, we have Af 



Npla 6 ; here N = N A - 



-Nb, Po = ^APo,A + (1 - 4>a)po,b, 
I, / = N A /N, a A (a s ) denotes 



a = y/fa 2 A + (1 - f)a 

the statistical segment length of the A (B) block, po A 
(/°o,b) the A (B) monomer number density, and </>a the 
overall volume fraction of the A block. Fig. [2] shows 
Af vs. the bulk lamellar period Lq for nearly sym- 
metric DBC melts commonly used in experiments 17]; 
the range of Lq = 10 ~ lOOnm, which is of interest 
for most applications of DBC, roughly corresponds to 
Af = 500 ~ 20,0000. 

In all conventional molecular simulations with hard 
excluded- volume interactions (e.g., the Lennard- Jones 
potential or the self- and mutual-avoiding walks), AT is 
on the same order of magnitude as N, which is at most 
one hundred or so for concentrated polymer solutions 
or melts (e.g., the leftmost data point in Fig. [JJ is at 
Af ~ 139). Using soft potentials that allow particle over- 
lapping is therefore the only way to study DBC melts 
within the above A^-range (i.e., with experimentally ac- 
cessible fluctuations) at present, where nN — pqV no 
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FIG. 2: (Color online) The invariant degree of polymerization 
Af vs. the bulk lamellar period Lo for nearly symmetric DBC 
melts commonly used in experiments. [l7l. [lq| 



longer holds and N becomes a chain discretization pa- 
rameter that does not correspond to the actual chain 
length used in experiments. This point is crucial for un- 
derstanding coarse-grained models with soft potentials, 
the use of which is the basic idea of the recently pro- 
posed fast MC simulations [ToL [l9j . 

In this work, we perform extensive FOMC simulations 
to quantify the ODT shift of symmetric DBC melts 
with the following Hamiltonian: H = W c + H E , 

where U c = (3k B T/2a 2 ) £Li (Rfc.^f-i ~ ^k, s f 

is the Hamiltonian due to chain connectivity, and 
H E (1/2/spo) / drdr' [p A (r) + p B (r)] u„(|r 

r'|)[p A (r') + /5 B (r')] + (x/po) J drdr'p A (r> (|r - 
r'|)/9B(r') is the Hamiltonian due to non-bonded in- 
teractions; here a is the effective bond length, Rfe jS 
denotes the spatial position of the s th segment on 
the k th chain, k and x are the generalized Helfand 
compressibility 20] for the copolymer melts and the 
generalized Flory-Huggins interaction parameter for 
the repulsion between A and B segments, respectively, 
Po = nN/V, p A (r) " ^" 
PbW = ELiEf 

number density of A and B segments at spatial position 
r, respectively, and u (r) = (15kBT/2na 3 )(l — r/cr) 2 
for r < a and otherwise is a normalized (i.e., 
J druo(|r|) = fcgT) isotropic pair potential depending 
only on the distance r between two segments with a 
denoting the finite interaction range. Note that Uo(r) is 
essentially the same as the potential for the conservative 
force used in dissipative particle dynamics [2_l|; its 
Fourier transform is positive definite, thus avoiding the 
formation of clustered crystals at large po[12|- Taking 
i? e ,o = \/N — la as the length scale, we have five 
parameters in our model: Af, xN, N/k, N, and a/a. 

Our simulations are performed in a canonical ensemble 
with trial moves of random hopping, reptation, pivot [23j. 
and box-length change. Note that the highly efficient 
pivot algorithm cannot be used in multi-chain Simula- 



= ELi£;=V(r-R*, s ) and 

=n a +i ^( r — R*,s) are the microscopic 
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tions with hard excluded-volume interactions due to its 
extremely small acceptance rates. With soft potentials, 
however, we can achieve an acceptance rate of about 70% 
here. We also use replica-exchange [lij at different xJV 
(with acceptance rates of 50^80%) to further improve 
our sampling efficiency. 

In simulations of periodic structures such as lamellae, 
the periodic boundary conditions limit the allowed ori- 
entations of the structure and thus its period. In par- 
ticular, for lamellae with a normal direction n in a sim- 
ulation box with length Lj in the j(= x, y,z) direction, 
Lj] • n = njL(n) must be satisfied, where nj is the num- 
ber of periods contained in the box along the j direc- 
tion (which could be 0) and L the lamellar period; this 

gives L(n) = l/ (m/L,-) 2 .[25| For a fixed-length 
box as commonly used in canonical-ensemble simulations, 
both the lamellar orientation and its period are therefore 
discretized (i.e., L can hardly be L ); this problem is 
the most severe for cubic boxes. |25l] To eliminate it, we 
change box lengths at a fixed V (i.e., using rectangular 
boxes) with segmental positions re-scaled according to 
the Metropolis acceptance criterion. Compared to sim- 
ulations in an isothermal-isobaric ensemble, this has the 
advantage that N is constant during each simulation. 

We use a new order parameter for ODT, which char- 
acterizes the degree of positional order in lamellae. For a 
given direction n, we calculate the volume fraction of 
A segments <J>a(£) = PA.(t)/\pA.{t) + pB (t)] as a func- 
tion of position t along n in a collected configura- 
tion (averaged over directions perpendicular to n), and 
V>(n) = |J dt exploit/ L(n)]f(t)\ / J df/(f) with f(t) = 
1 — \2<pj^(t) — 1| and i = \f—l. The order parameter 

€ [0,1] for a collected configuration is defined as the 
largest ip over all possible lamellar orientations in the 
box. Fig. [3] shows the ensemble average of <!», (ty), and 
d (vfr) /d(yJV) vs. xN. We further use multiple histogram 
reweighting[26| to accurately locate ODT according to 
the equal-weight criterion[27| described below. 

The inset of Fig. [3] shows the histogram of f at a 
given x, P 6 (*;x), with d^P b (^;x) = 1- The double 
peak near ODT is a signature of the first-order phase 
transition. The ODT is determined as x* at which 
d*Pb(*;x*) = / ( «>d¥A(*;X*)- To estimate the 
statistical error of x* > we calculate ODT using the first- 
and second-half of our samples collected after equilibra- 
tion, respectively, and take three times their largest de- 
viation from x* determined using all the samples as the 
error bar. Note that, even with the same order param- 
eter, there are other ways to determine the transition 
point [Hj], which give slightly different ODT for finite sys- 
tems as shown in Fig. [3] Systematic study of the finite- 
size effects requires much more extensive simulations and 
is out of the scope here. 

With the same Hamiltonian as used in our FOMC sim- 
ulations, we determine the mean-field ODT Xmf usm g 
the random-phase approximation TableUlists Xmf-^ 
for various JV and a /a. We note that the mean- field re- 
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FIG. 3: (Color online) Ensemble average of the order param- 
eter, (*}, and d(*) /d(xN). The error bar of (*) is taken 
as three times the standard deviation with sample correla- 
tions taken into account. The vertical line marks the ODT 
X*N determined using the equal- weight criterion[23], which 
is very close to the location of the maximum of d (9) /d(xN) 
(i.e., another way to determine ODT[2£|). The inset shows 
the reweighted histogram of ^ at x*> Pb(^',X*)i where the 
vertical line marks (*) (x*). N = 10, a/a = 0.3, N/k = 50, 
and JV = 10, 000. 



TABLE I: Mean-field ODT Xmf^- 







a/a = 0.3 


a/a = 0.lVl§ 


a/a = 2/V3 


N = 


10 


10.047 




11.427 


N = 


20 


10.405 


10.462 


11.102 



suits are independent of N/k and M, and can be under- 
stood based on our previous workfiol [29|. 

Collecting all the data, Fig.2]shows our results of ODT 
shift vs. N at various JV", a/a and Nf k. [30j We find 
X*/xmf ~ 1 fx A/" _fc in all the cases. While this is con- 
sistent with the functional form of FH prediction, our 
ODT shift is larger than their prediction for all the cases 
we have studied. We also find a decrease in the nega- 
tive exponent k around A7 cp = 2(JV - l) 3 /N 2 (a/a) 6 in 

the cases of a/a — 0.3 and 0.1-\/T9 ~ 0.44, which corre- 
sponds to the FCC close packing of polymer segments as 
hard spheres. This is therefore due to the local packing 
of segments (i.e., the short-range correlations). Denoting 
the fc-value for Af < A7 cp by k\ and that for Af > A7 cp 
by k 2 , we find that ki — k 2 decreases with increasing N 
(i.e., in the case of a/a — 0.3), which leads to decreasing 
a (thus weaker correlations) at constant i? e ,o- One may 
therefore expect a single power- law decay of X* /Xmf ~ 1 
with increasing M in the limit of JV — > oo (i.e., a —> 0), 
which has no such correlation and is the case studied by 
Fredrickson and Helfand 14] . On the other hand, increas- 
ing a /a at constant JV leads to stronger correlations and 
thus increasing k\ — k 2 (i.e., in the case of JV = 20). We 
also note that varying JV at constant a /a exhibits both 
the correlation and chain discretization effects. To exam- 
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FIG. 4: (Color online) Our results on the ODT shift. The 
negative slope k is given next to each line. See text for details. 

ine the latter alone, one may compare the case of N = 10 
and a/a = 0.3 with that of N = 20 and a/a = 0.1y/l8 
(both of which have a = Q.l-R e ,o)- We expect that the 
chain discretization effects diminish quickly with increas- 
ing N, as supported by Fig. 6(a) of Ref. [1(| showing how 
XmfN varies with N. 

In the case of a/a — 2/Vo s» 1.15, Af cp < 15 for both 
iV-values, and we did not perform simulations at N < 



100. This large a /a- value also makes simulations at M > 
2000 expensive to do. Nevertheless, we see from Fig. 0] 
that increasing N increases the ODT shift in this case, 
in contrast to that of a/a = 0.3. On the other hand, as 
we increase N/k from to 50, /c2 decreases here and the 
two lines with the same N in Fig. |4] cross at Af CI ~ 1300 
(for N = 10) or 2500 (for N = 20). Therefore, x*N 
decreases for Af < Af CI and increases for N > J\f CT with 
increasing N/k. While the latter is consistent with the 
case of a/a = 0.3, we find that both k\ and k-i increase 
with increasing N/k there. 

To summarize, we have performed extensive FOMC 
simulations with experimentally accessible fluctuations 
to systematically and unambiguously quantify the ODT 
shift of symmetric DBC from the mean-field prediction. 
Exactly the same model system (Hamiltonian) is used 
in both our simulations and mean-field theory; the ODT 
shift is therefore due to the fluctuations/correlations ne- 
glected by the latter. We hope this work will stimulate 
the development of advanced theories better describing 
the fluctuation/correlation effects identified here. 
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work was supported by the U.S. Department of En- 
ergy, Office of Basic Energy Sciences, Division of Mate- 
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